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ABSTRACT 

We discuss the skeleton as a probe of the filamentary structures of a 2D random field. 
It can be defined for a smooth field as the ensemble of pairs of field lines departing 
from saddle points, initially aligned with the major axis of local curvature and con- 
necting them to local maxima. This definition is thus non local and makes analytical 
predictions difficult, so we propose a local approximation: the local skeleton is given 
by the set of points where the gradient is aligned with the local curvature major axis 
and where the second component of the local curvature is negative. 

We perform a statistical analysis of the length of the total local skeleton, chosen 
for simplicity as the set of all points of space where the gradient is either parallel 
or orthogonal to the main curvature axis. In all our numerical experiments, which 
include Gaussian and various non Gaussian realizations such as fields and Zel'dovich 
maps, the differential length is found within a normalization factor to be very close 
to the probability distribution function of the smoothed field. This is in fact explicitly 
demonstrated in the Gaussian case. 

This result might be discouraging for using the skeleton as a probe of non Gausian- 
nity, but our analyses assume that the total length of the skeleton is a free, adjustable 
parameter. This total length could in fact be used to constrain cosmological models, 
in CMB maps but also in 3D galaxy catalogs, where it estimates the total length of 
filaments in the Universe. Making the link with other works, we also show how the 
skeleton can be used to study the dynamics of large scale structure. 



■ 1 INTRODUCTION 

The observed large scale distribution of galaxies presents re- 
markable structures, such as clusters of galaxies, filaments, 

■ sheets and large voids. It is widely admitted that these struc- 
' tures grew from small initial fluctuations through gravita- 
, tional instability. At very large scale, the filamentary pattern 

seen in the cosmic web is expected to be similar to that of 
the initial field (e.g., Bond, Kofman & Pogosyan 1996). Since 
these primordial inhomogeneities also imprinted the temper- 
ature fluctuations seen now in the Cosmic Microwave Back- 
ground (CMB), the characterization of the observed large 
scale structures both in galaxy catalogs and in CMB maps 
can help to probe the nature of these primordial fluctua- 
tions, in particular whether they are Gaussianly distributed 
or not. Furthermore, a rigorous topological description of 
the observed structures is necessary to constrain efliciently 
models of large scale structure. For instance, a precise defl- 
nition is needed for clusters of galaxies before inferring any 
constraints from their studies, e.g density and temperature 
proflles but also luminosity function and clustering. Analo- 
gously, a precise and practical deflnition of fllaments would 
allow us to use them similarly as cluster of galaxies. 

Various methods have been proposed to characterize 
the morphology of large scale structures. In general, one 



studies the topological properties of excursion sets, i.e. re- 
gions below or above a density threshold. The statistics 
most explored up to know are the genus or the closely 
related Euler characteristic (see, e.g., Doroshkevich 1970; 
Gott, Melott & Dickinson 1986), the more complete set 
of observables given by Minkowski functionals (see, e.g., 
Mecke, Buchert & Wagner 1994) and related statistics such 
as shape finders (e.g., Sahni, Sathyaprakash & Shandarin 
1998), but also estimators based on percolation analysis 
(see, e.g., Zel'dovich 1981; Zel'dovich, Einasto & Shandarin 
1982; Shandarin 1983; Bhavsar & Barrow 1983) and mini- 
mum spanning tree construction (e.g., Barrow, Bhavsar & 
Sonoda 1985), such as for example structure functions or 
related shape estimators based on moments of inertia (see, 
e.g., Babul & Starkman 1992). 

In general, topological descriptors are primarily used to 
constrain the level of non Gaussianity in the sample, since 
there often exists analytical predictions in the Gaussian case 
(see, e.g., Doroshkevich 1970 for the genus; Tomita 1986 for 
the Minkowski functionals). To do so, other estimators exist 
also, based on spatial correlation analysis, such as higher- 
order correlation functions in real or Fourier /harmonic space 
(e.g., Peebles 1980 and references therein), the probability 
distribution function (pdf) of the smoothed density fleld and 
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its moments (e.g., Bernardeau et al. 2002 and references 
therein), higher order moments the wavelet coefficients (e.g., 
Aghanim & Forni 1999; Hobson, Jones & Lasenby 1999), 
peak and excursion set subcomponents statistics (see, e.g., 
Bardeen et al. 1986; Bond k Efstathiou 1987 and Dore et 
al. 2003 for a somewhat related statistic), phase correlation 
analysis (e.g., Chian & Coles 2000; Naselsky, Novikov & Silk 
2002), etc. In principle, all these statistics combine the data 
in very specific ways, so they altogether provide complemen- 
tary analysis of the data.^ However, at variance with tradi- 
tional statistical estimators, topological estimators help as 
well to characterize in a very intuitive way the topology of 
structures in terms of filaments, sheets, clusters and voids 
quantitatively. For instance, the Minkowski functionals pro- 
vide a complete basis of simple estimators to estimate the 
morphology of an excursion set (see, e.g., Kerscher 2000 for 
a review on the subject), e.g. its degree of compacity and of 
filament arity. 

In this paper we focus on the skeleton, which aims at ex- 
tracting from the cosmic web its filamentary pattern. More 
specifically, the goal is to draw in the observed structure a 
set of lines which reproduces well the filamentary pattern 
guessed by eyes. A natural tool to do so in set of points such 
as galaxy catalogs is the minimum spanning tree (e.g.. Bar- 
row, Bhavsar & Sonoda 1985). It is a connected structure 
superposed to the set of points, with no loop and which is 
the shortest possible. Of course, using the minimum span- 
ning tree as such is not very helpful since it is highly irregular 
and it is difficult to establish a link to the large scale fea- 
tures of interest, but there are technics to filter small scale 
noise consisting in "pealing" the tree, i.e. removing from 
it short branches. Even if it is successful in extracting the 
main filamentary features from the catalog,^ the minimum 
spanning tree remains by nature unsmooth and difficult if 
not impossible to manipulate in order to perform analytic 
calculations. 

The technique we aim to employ in this paper to extract 
the skeleton from the data sample is completely different 
and relies on Morse theory (see, e.g., Milnor 1963; Colombi, 
Pogosyan & Souradeep 2000; Jost 2002). It requires the field 
to be sufficiently differentiable and non degenerate as ex- 
plained more in details later and thus some smoothing with 
e.g. a Gaussian window of the data file is necessary to use 
such technique. "i" ^ 

^ see, e.g., Shandarin 2002 for a comparison of the pdf and the 
Minkowski functionals as estimators of non Gaussianity in 2D 
maps; Phillips & Kogut 2001 for a comparison between genus, 
extrema correlation function and bispectrum; Barreiro, Martmez- 
Gonzalez & Sanz 2001 for a comparison between number density, 
eccentricity and Gaussian curvature of hot spots, and genus as 
estimators to probe non-Gaussianity in CMB samples, 
t see, e.g. Fig. 1 of Doroskevich et al. 2001 for a nice illustration 
on the Las Campanas Redshift Survey. 

t note that the concept of smoothing introduce a scale in the 
problem: smoothing at different scales will not produce the same 
skeleton but will have interesting links to dynamics as discussed 
in the end of this paper. 

§ Hence, at variance with the minimum spanning tree method, 
which has the advantage to deal directly with the discrete nature 
of galaxy sample, our method will be difficult to apply to real 
galaxy catalogs unless smoothing at sufficiently large scales. 



We shall see that the skeleton can then be then rigor- 
ously defined as a set of pairs of special field lines departing 
from saddle points. The problem is that the skeleton defined 
as such is non local: indeed, to draw any field line, one has 
to resolve the trajectory of a particle following the equation 
of motion given by dr/dt = Vp. This non local nature of 
the skeleton makes analytic predictions rather difficult. Fur- 
thermore, as discussed in Appendix A, it is difficult to find 
a reliable algorithm to draw it. 

The main points of this paper, which focusses on the 
2D case, are the following: (i) find a local approximation of 
the skeleton to address the issues just raised above, (ii) test 
this local approximation as a statistical tool to probe non 
Gaussian features of the density field in e.g. CMB maps, (iii) 
establish the link to dynamics. The last point will be only 
treated superficially through simple illustrative examples, 
relying mostly on the Zel'dovich approximation (Zel'dovich 
1970), since it is clearly more interesting to treat it in detail 
in the 3D case. 

This paper is thus organized as follows. In Section 2, 
We define the skeleton in the framework of Morse theory 
and discuss some of its properties. We find a local approxi- 
mation to it, relying on two independent methods. This local 
skeleton is contained in the set of points of space where the 
gradient of the density field is aligned with one axis of the lo- 
cal curvature, that we call the total local skeleton. We show 
through examples that the local skeleton and the real skele- 
ton are quite alike, both by visual inspection and by com- 
paring their respective lengths and we discuss the differences 
found. In Section 3, we study the differential length of the 
total local skeleton as a function of density threshold. We 
find experimentally that it scales very much like the pdf of 
the smoothed density field, a property that we demonstrate 
explicitly in the Gaussian case. Finally, Section 4 discusses 
the results and makes the link to dynamics. An extensive 
appendix discusses the numerical calculations of this paper, 
which were performed with a dedicated package. 



2 THE SKELETON OF A 2D RANDOM FIELD 

In a two-dimensional field, one would naturally define the 
skeleton as a set of ridges connecting local maxima and sep- 
arating under-dense regions. In what follows, we first give 
the practical mathematical definition corresponding to this 
view (§ 2.1). It is shown that the skeleton is a set of special 
field lines, i.e. a particular set of curves parallel to the gra- 
dient of the field, and passing through maxima and saddle 
points. However this definition is not very practical at least 
from the theoretical point of view, because it is non local 
and makes analytic predictions difficult. To enforce locality, 
we define an other set of curves that is aimed to be close to 
the real skeleton (§ 2.2). To do so we use two different ap- 
proaches, which actually end in the same definition for the 
local skeleton. The first one consists in Taylor expanding the 



' Note however, that dynamics in 2D can be still of interest, 
particularly in relation to reconstruction of the projected mass 
distribution in weak lensing experiments: for instance the skeleton 
can be used as a tool to test the quality of the reconstruction 
methods. 
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Figure 1. An example of Gaussian field and its smoothed counterpart. 

Left panel: a periodic realization of a 2D Gaussian random field on a grid of size 1024 x 1024 pixels with a scale-free power-spectrum, 
verifying P{k) oc k'^ with n = —1. 

Right panel: the same field smoothed with a Gaussian window of radius 25 pixels. 



special field lines in the neighborhood of saddle points and 
local maxima while the second one consists in finding points 
along isocontour lines such that the gradient of the density 
field has extremal magnitude. In § 2.3, our arguments will 
be illustrated by practical examples on a Gaussian field and 
its Zel'dovich mapping. As a probe of the local skeleton, in 
addition to visual inspection, we shall compare, for the ex- 
amples considered here, its length as a function of density 
threshold with the length of the real skeleton. 

In what follows, we consider a 2D random field, p(r), 
e.g. a temperature map of the CMB. To assume sufficient 
differentiability we convolve it with a gaussian window of 
size i, as illustrated by Fig. 1. The smoothed field, still noted 
p, is furthermore supposed to be sufficiently non degenerate, 
and in particular has the following properties: 

(i) Its gradient cancels in a discrete set of critical points, 
which can be separated into three subclasses, local minima, 
local maxima and saddle points; 

(ii) The eigenvalues Ai > A2 of its Hessian, ?i = 
d'^p/dridrj, verify the following properties: the regions of 
space where Ai = or where A2 = are sets of smooth 
curves never passing through critical points. The intersec- 
tion of these two sets of curves, where Ai = A2 = 0, is 
therefore a discrete set of points not containing any critical 
point. 

Given these last definitions for Ai and A2, the local maxima, 
saddle point and local minima verify respectively > Ai > 
A2, Ai > > A2 and Ai > A2 > 0. 

2.1 Definition 

We first define the peak patches (void patches) as the regions 
of space containing all the points converging to the same 
local maximum (local minimum) while going along the field 
lines in the direction (opposite direction) of the gradient, 
Vp = dp/dvi. 

The skeleton (of over-dense regions) is defined as the 
borders of the void patches (and a dual skeleton can be 
similarly defined as the borders of the peak patches). It is 



easy to show that it passes through all the saddle points and 
the local maxima. It would be out of the scope of this paper 
to go further in the mathematical details of the topology 
of the skeleton, but one can list the following well known 
properties, valid only if there is no unexpected degeneracies 
(e.g., Jost 2002), and which can be easily verified by visual 
inspection of Fig. 2 (top left panel): 

• The nodes of the skeleton are the local maxima, where 
multiple lines of the skeleton can converge. In general, be- 
cause Ai > A2, these lines tend to converge along the axis 
aligned with the eigenvector associated to Ai (Fig. 3, left 
panel) . 

• Two local maxima cannot be directly connected to- 
gether, there is always a saddle point in between; 

• Saddle points cannot be nodes of the skeleton. Indeed, 
there are only four field lines connected to each saddle point: 
two unstable fields lines arriving from opposite directions, 
locally parallel to the eigenvector corresponding to A2 < 0, 
and two stable fields lines departing in opposite directions 
locally parallel to the eigenvector corresponding to Ai > 
(Fig. 3, right panel). These too last field lines locally coincide 
with the skeleton and end to a local maximum. 

From the last argument, the skeleton can be seen as the 
ensemble of pairs of stable fields lines departing from saddle 
points and connecting them to local maxima. The skeleton 
field lines can thus be drawn by going along the trajectory 
with the following motion equation 



starting from the saddle points, and with initial velocity par- 
allel to the major axis of the local curvature (i.e. parallel to 
the eigenvector of the Hessian corresponding to Ai). Tra- 
jectory is followed until convergence to a local maximum. 
This procedure was actually used to draw the skeleton, as 
explained in details in Appendix. 
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Figure 2. Skeleton and its local approximation for the Gaussian field of Fig. 1. 

Upper left panel: the skeleton is drawn as well as the critical points: local minima in yellow, saddle points in orange and local maxima 
in red. As discussed in the text, the skeleton passes through all the maxima and the saddle points. The local maxima are the nodes 
where several lines converge, while the saddles points have only one line passing through. Note as well that local maxima are always 
connected to saddles and reciprocally, except in e.g. the lower left of the panel, when one can see three saddle connected to each 
other. This configuration is theoretically forbidden unless there is some degeneracy in the field, that we suspect is due to our numerical 
implementation, as further discussed in Appendix. 
Upper right panel: the skeleton is superposed to the smoothed field. 

Middle left panel: same as for the upper left panel, but for the local approximation of the skeleton. The dark plus light blue lines assume 
S = [eq. (8)], while the light blue lines verify the more constraining conditions given by eqs. (2) and (3). 
Middle right panel: same as upper right panel but for the local approximation of the skeleton. 

Lower left and lower right panels: the local approximation and the real skeleton are again superposed to the smooth field, but restricted 
to over-dense regions p > (p). 
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Figure 3. Topology in the neighborhood of critical points (inspired from Jost 2002). 

Left panel: expected topology of field lines nearby a maximum, if Ai > A2. Except for two vertical field lines along the eigenvector 
corresponding to the smallest eigenvalue, A2, all the lines converging to the node tend to be aligned with the horizontal axis, corresponding 
to Ai. If Ai = A2, we would face a degenerate situation where all the directions are equivalent. 

Right panel: expected topology of field lines nearby a saddle point, if Ai > > A2. There are only four field lines connected to this 
point, aligned with the two eigenvectors. The two vertical field lines, corresponding to A2, are unstable since all the other field lines are 
diverging away from them. Reversely, the two horizontal lines corresponding to Ai are stable lines. 



2.2 Local approximation 

As discussed in Appendix, equation of motion (1) is not 
easy to solve in practice, even for a smooth field sampled on 
a finite but thin grid. Furthermore, analytic predictions are 
very difficult since eq. (1) is non local. This motivates the 
need for an approximation of the real skeleton with a local 
criterion on the density field and its derivatives of various 
orders. We shall do so by two means, the first using a more 
mathematical approach, the second using a more physical 
approach. 

The mathematically motivated derivation consists in 
Taylor expanding the field around its saddle points and 
its local maxima. On the skeleton nearby these points, we 
clearly have, at leading order (and except for degenerate 
cases) 



A2 < 0, 



(2) 
(3) 



This entails a natural definition for the local skeleton: it 
consists of any point of space where eqs. (2) and (3) are 
verified. 

The physically motivated approach consists in consid- 
ering the field as a landscape, where the third coordinate, 
rs, is given by rs — p(ri^r2)- In that case, the isocontour 
lines of the density field are natural objects for the analyses. 
Let us consider two pieces of isocontour lines A and B very 
close to each other and let us move from A to following 
the gradient. We expect the skeleton to take either the short 
possible or the longest possible path between A and B. Since 
the path length is inversely proportional to the magnitude 
of the gradient, the points of interest are those where |Vp| 

is locally an extremum along the isocontour.il 



An accurate examination of the neighborhood of local max- 



This translates mathematically as follows. If we denote 
s a curvilinear coordinate along an isocontour, (ri(s), r2(s)), 
we have by definition. 



dp dri ^ dp dr2 
dri ds dr2 ds 



0, 



(4) 



with the normalization 
2 / X 2 



The gradient is locally an extremum along the isocontour if 



(IVpH^o, 



dt 

which gives, using eq. (4) 
S = 



(6) 



dp_dp_ f d^p d^p 

dri dr2 \ dri 

d^p (\dp-^ 
dridr2 



dr2 



dp 



dri 



= 0. 



In fact, it is fairly easy to rewrite this equation as 
S = det (HVp, Vp)=0. 



(7) 



(8) 



Therefore, the condition <S = is equivalent to say that the 
gradient is an eigenvector of the Hessian. 

However, there is a supplementary condition which 



ima and saddle points suggests that the skeleton should take the 
longest possible way, which in fact implies that the magnitude of 
the gradient is a local minimum along the contour line. However, 
enforcing such a condition would lead us to examine expressions 
involving third derivatives of the density field, in disagreement 
with a leading order approach. Instead, we are going to use less 
realistic but simpler criteria on the local curvature of the isocon- 
tour lines nearby the extrema of the density gradient magnitude 
to select the points of interest. 
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comes out naturally: while walking from one field line to 
another, one prefers to stay on a ridge, that is on the points 
where the curvature of the isocontour is positive, i.e. 



^ - |Vp| ds^ ^ 



(9) 



which translates in, after some algebra based on eqs. (4) and 
(5), 



C = 



|Vp|3 



\7p±n\7p± > 0, 



(10) 



where Vpj_ = {dp/dr2, —dp/dri). So, we have to select 
among the points verifying eq. (8) those which have C > 0. 
Since, Vp is an eigenvector of H, so is Vp± , therefore 



A2 
■|Vp| 



Ai 

IVpl' 



(11) 



After a simple examination of the various cases, Ai > A2 > 
0, Ai > > A2 and > Ai > A2, we finally obviously 
converge again to eqs. (2) and (3), except for hills, > Ai > 
A2, where equation (10) allows the gradient to be aligned 
with both axes of the curvature. In this last case, we can see 
however that the situation TiVp = A2Vy9 contradicts the 
"natural" definition of a ridge. Such a ridge would indeed 
be more curved along its path than orthogonally to it, a 
situation clearly unrealistic in the neighborhood of a local 
maximum. 



2.3 Examples 

We now compare the local and real skeleton by visual in- 
spection of Figs. 2 and 4, which respectively correspond to 
a Gaussian realization and its Zel'dovich mapping. We also 
measure the length of the skeletons as a function of threshold 
for these two particular examples. 

Clearly the local skeleton is an excellent approximation 
of the real one. Indeed most of the large scale features are 
very well captured, particularly in the vicinity of maxima 
and saddle points, as a result of our perturbative approach. 
In agreement with intuition, the more filamentary is the 
field, the better is the agreement: the local skeleton seems to 
perform better in the Zel'dovich map than in the Gaussian 
one. However, connectivity of the local skeleton is not en- 
sured, at variance with the real one. Furthermore, there are 
little spurious structures in void patches that do not match 
any line of the real skeleton. By restricting the comparison 
to over-dense regions, a large number of these structure dis- 
appear and the agreement improves significantly, at least 
visually. 

A more quantitative analysis can be conducted by com- 
paring the measured length of the real and local skeletons in 
regions where the density exceeds a given threshold, as illus- 
trated by Fig. 5. Contrary to what would suggest the visual 
inspection of Figs. 2 and 4, the local skeleton is systemat- 
ically slightly shorter than the real one. The total lengths 
differ by about 20 percents, both for the Gaussian smooth 
field and its Zel'dovich mapping. However, as illustrated by 
left panel of Fig. 3, skeleton field lines converging to a lo- 
cal maximum tend to superpose along the major axis of the 
local curvature, which produces multiple lines. This feature 
inherent to the Lagrangian nature of the real skeleton is 



missing in the local skeleton, due to its local, Eulerian na- 
ture (see also Appendix). Hence, it is not surprising that the 
local skeleton is shorter than the real one. 



3 LINK TO STATISTICS 

In this section we focus on the local skeleton. For the sake 
of simplicity, we study from now on the total local skeleton, 
defined as the full set of points satisfying the condition S = 
[eq. (8)]. We first examine the Gaussian case in § 3.1, where 
specific analytic results are derived and then confronted to 
numerical experiments. The normalized differential length 
of the skeleton as a function of density threshold is seen to 
be very close to the probability distribution function (pdf) 
of the smoothed field, i.e. a Gaussian. 

We thus consider in section § 3.2 some examples of non- 
Gaussian fields, namely distributions with n degrees of 
freedom, the Zel'dovich mapping discussed previously and, 
finally, an extreme case where the density contrast is lo- 
cally enhanced along lines with random orientations. In all 
these cases, we find again that the differential length of the 
skeleton scales very much like the pdf, with slightly worse 
noise properties as expected since this estimator relies on 
derivatives of the field. This intriguing result, despite its 
mathematical beauty, might look discouraging for using the 
skeleton as a test of non-gaussianity. However, our analysis 
does not use the supplementary information provided by the 
skeleton, namely its total length, which is considered here 
as an arbitrary normalization. 



3.1 The Gaussian case 

In this section, after defining a set of useful notations 
(§3.1.1), we compute a general expression for the differential 
length of the skeleton, /, as a function of density threshold 
(§3.1.2). Then we concentrate on the Gaussian case (§3.1.3) 
and show that the shape of / depends only on a single spec- 
tral parameter, 7, defined below. We demonstrate that for 
7 = 0, / is exactly given by a Gaussian. Given the level of 
complexity of the calculations, we are however lead to rely 
on a power expansion in 7 around a Gaussian to examine 
the case 7 > 0. This power expansion is checked carefully 
against numerical experiments (§3.1.4), which indeed show 
that / deviates only weakly from a Gaussian. 



3.1.1 Spectral parameters and dimensionless variables 

From now on, we assume without loss of generality that the 
random smoothed field, p has zero average {p) — 0. For 
convenience we define the following spectral parameters: 



2 
2 



2(p?)=2(pi), 



0-2 = 3(^11) = 3(^22) = 8(P22), 



(12) 
(13) 

(14) 



y (Piii) = y (P222) = 16(pfi2) = 16(pf22), (15) 



7 = 0-1/(0-00-2), 

7 = 0-2/(o-lO-3), 



(16) 
(17) 
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Figure 4. Skeleton and its local approximation for the Zel'dovich mapping of the smoothed Gaussian field of Fig. 1. 

Upper left panel: Zel'dovich mapping of the smoothed field: the Lagrangian displacement field P was normalized so that V.P = —(p — 

{p))/a{p), where is the variance of the initial smooth map. 

Upper right panel: the field of the top left panel smoothed with a Gaussian window of radius 12.5 pixels. Such smoothing is necessary 
to get rid of caustics and to enforce the differentiability required to measure the skeleton. The smoothing scale is such that the large 
scale features outside caustics are preserved: it has to be small compared to the initial smoothing radius of 25 pixels and large enough 
compared to the pixel size. 

Middle left and middle right panels: respectively, the real and the local skeleton superposed to the Zel'dovich map. 

Lower left and lower right panels: respectively, the real and the local skeleton superposed to the Zel'dovich map, but restricted to 
over-dense regions. 
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Figure 5. Real versus local skeleton length comparison. 

Left panel: the skeleton length in units of sample box length as a function of density threshold, measured for the real (thin curve) and the 
local (thick curve) skeleton in the smooth Gaussian field of Fig. 1. S/a = (p — {p))/a(p) is the density contrast in units of the variance 
of the smoothed field. 

Right panel: same as left panel but for the smoothed Zel'dovich map (top right panel of Fig. 4). 



where pi = dp/dn, pij = d'^p/dvidrj and pijk = 
d^p/dndrjdrk {ij^k = 1,2) are its gradient, hessian and 
matrix of third derivatives, respectively. Using these param- 
eters one can consider the 10 following dimensionless vari- 
ables: 

X = , Xi = -, Xij = — , Xij}^ = — , (1^) 
(To (Ti (72 as 

and the dimensionless function s = S/{cf\cf2)- That is, ac- 
cording to eq. (7): 

S = X1X2(X11 - X22) + Xl2{xl - xl). (19) 

Thus, the points of the random field where the first and 
second derivatives satisfy the condition s = define the 
total local skeleton. Finally, the derivatives 



Si = Si/{(7i(7s) 



(20) 



will also be useful. It is worth noting here that Si depends 
on 7, but not on 7, i.e. si = Si{xi,Xij,Xijk,'y). 



3.1.2 Length of the skeleton: general expressions 

If we denote Vs{p,S, 81,82) the joint probability distribu- 
tion of the variables p, <S, and <Si for i = 1, 2, the expected 
average of the skeleton length C(pth) per unit area above 
some threshold pth can be derived the following way. 

Let us consider a straight line along the direction n and 
which intersects the isocontour lines, 8 = 0, at some point, 
where p > pth- In the vicinity of such a point, where 8 = 
and d8 = 8idri, we can integrate Vsdpd8d8id82 over dn 
from —dri/2 to dri/2 and we get, obviously: 

J d8d8id82dpVs{p, 0, 81,82) = 

rie[-dri/2,dri/2], p>Pth 

dn J \8i\d8id82dpVs{p,0, 81,82). (21) 
p>Pth 

This integral represents the probability to find the point 
8 = along the line r2 =const in the range [n — dri/2, n + 



dr 1/2]. Note that the absolute value of 81 is here necessary 
since we want to take into account both up-crossing and 
down-crossing points. Then, the elementary length of the 
isocontour line 8 = inside the square [ri — dri/2,ri -h 
dri/2]r2 — dr2/2,r2 + dr2/2] is dr2/ cos(q;), where a is the 
angle between r2 and the isocontour fine. Since cos(q;) = 
\8i\/ ^ 81 + <S|, using eq. (21) one gets 

C{pth)dridr2 = 

dridr2 J (ip6Z5icZ52V^Sf+5| P.(p,0,<Si,<S2). (22) 

p>pth 

This equation represents the average length of the skeleton 
per element of area dridr2. In terms of dimensionless vari- 
ables, it rewrites 

C{xth) = / dxdsids2 — v/sj -h s| Ps(a:,0, si, S2). (23) 
J ^2 

We shall now look for an analytical expression for this func- 
tion C in the case of a Gaussian field. 



3.1.3 Towards an analytic expression for a Gaussian field 

Deriving eq. (23) we did not consider any special features 
of the joint probability function Vs^ thus this equation is 
true for any random field. Unfortunately the derivation of 
function Vs is not easy, even in the Gaussian case, that we 
examine now. 

We consider the 10 components random vector, a. 



a. — {x, Xi, Xij, Xijk) {i^ j-) k — 1,2), 



(24) 



and the probability distribution function P(a). In the Gaus- 
sian case, this can be written as 



P(a) 



[(27r)io|M|]-5 



exp 



--aM-^a^ 



(25) 



where M is the covariance matrix, M = (aa^), and |M| its 
determinant. Then, eq. (23) can be rewritten using 7^ (a) 
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s|7'(a)da. 



(26) 



CC>CCth5 s = 



For further investigation, it is particularly convenient to con- 
sider the following variables: 



q = xii 



X22, 



U = 2xi2, V = X\i-\-X22- 



(27) 



We are interested in the distribution of the skeleton length 
over just one variable Xth- It is worth mentioning that using 
the variables defined in eq. (27), x correlates only with v 
and does not correlate with either g, xi nor xijk^ i.e. 



{xq) — (xu) = (xxi) = {xxijk) = 0, 
(xv) — —7. 



(28) 
(29) 



Taking into account eqs. (27) and (29), one can represent 
P(a) in the following way: 



^(a) 



: exp 



(x + 



2(1-7^) 



V(xi,q,u,v,Xijk)- (30) 



In order to take into account the condition s = in eq. (26), 
we should perform one more substitution: 



xi = r cos{(f), 
q = pcos{2^|;), 



X2 = rsm{(f), 



(31) 
(32) 



It is easy to see that with these new variables, the function 
s in eq. (19) reads: 



s = ^pr^ sin [2{(f — ip)'^ . 



(33) 



The condition s = is now transformed into a simple rela- 
tion between the angles (f and -0. The length of the skeleton 
can therefore be written as follows: 



27 



C(xth) = — dx dv ^ J=r / c/0 a/s? + 

<^2 J J ^/l^ J 



) (34) 

where c/Q = pdpr dr dxijk dip dip and is the usual delta 
function. This equation will give us the total length of the 
skeleton jCtot if we consider xth —00. From eqs. (19) and 
(34) one can see that the differential length normalized by 
the total length is: 



1 dC 

Ctot dx 



+00 

/ 



dv C{v) 



2(1-7^) 



(35) 



Therefore, remarkably, / is a function of x and only one 
spectral parameter 7. The quantity f{x,^)dx simply repre- 
sents the fraction of the skeleton length between the levels 
X and x-\-dx. However, the unknown function C{v) is rather 
cumbersome to estimate analytically. We examine in next 
section a way to avoid its calculation, but which relies partly 
on numerical experiments. 

3.1.4 Final expression in the case of a Gaussian field 

The first thing to notice, when examining eq. (35), is that 
in the limit 7 = 0, / is exactly a Gaussian: 

1 -.x2/9 



/(^,7 = 0) 



(36) 



We thus expect / to depart only weakly from a Gaussian if 
7 is small enough, which motivates for a power-expansion of 
/(x,7) in 7 around a Gaussian. 

According to eq. (35), /(x,7) satisfies the following par- 
tial differential equation: 

df _ d fdf 



^7 



7 ■ 



dx 



(37) 



Since /(x,7) is a probability distribution function, it should 
satisfy the following condition: 

Jf{x,^)dx = l. (38) 
To try to solve eq. (37), we examine solutions of the form 



/(^.7) = < ^^n(x)7^ 



27r 



Injecting this expression in eq. (38) leads to 
-h ngn = 0. 



d'^gn 



(39) 



(40) 



dx"^ dx 
Setting y = x/V2, we find 

d'^gn ^ dgn . ^ ^ .... 

_.2y—+2ng„ = 0, (41) 

a differential equation followed by Hermite polynomials, 
Hn{y). A more detailed examination of the possible solu- 
tions gives 

Cn + Dn I dy eMy^)/Hn{y) • (42) 



gn{x) = Hn{x) 



/X 



Questionable arguments based on enforcing the convergence 
of the moments of / with respect to x suggest Dn — 0. 
Symmetry of the Gaussian field with respect to x = implies 
C2n+i = 0. As a result we expect / to have the following 
form 



I n>0 J 



-x^/2 



27r 



(43) 



with Co = 1, from normalization (38). Note that this ex- 
pression is nothing but a Gram- Char lier expansion prior to 
standardization [see Stuart & Ord 1994, eq. (6.32)]. It can 
be valid in practice only if the departure from a Gaussian is 
weak. Alternatively, thus, we could have derived eq. (43) by 
trying to find solutions of the form 



/(a^,7) = I ^fe„(7)^2n(a;/V2) I -^e 



(44) 



using directly the Gram- Char lier expansion. 

The coefficients C2n, n > 1, remain to be determined. 
We therefore performed a set of numerical experiments 
which was used as well to test extensively our skeleton 
analysis software (see discussion in Appendix). We gener- 
ated scale-free random Gaussian fields with power-spectra 
P{k) oc /c"", n = 0, -1, and -2 (7 = 0.71, 0.58, 0.32, respec- 
tively). For each value of n we performed 100 realizations 
over a periodic grid of size 1024 x 1024 pixels (except for 
n = 0, where we used 2048 x 2048 pixels). The field was 
smoothed with a window of radius 5 pixels (10 pixels for 
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Figure 6. The measured differential length of the total local skeleton [eq. (8), left panels] and the pdf of the smoothed field for scale- free 
Gaussian random fields (right panels) with power-spectra P{k) k^, n = 0, —1 and —2 as indicated on each panel, as functions of 
normalized density contrast, rj = S/a = x = p/a. For each value of n, 100 realizations were performed on a periodic grid of size 1024 x 1024 
pixels (2048 x 2048 for n = 0), and then smoothed with a Gaussian window of radius 5 pixels (10 pixels for n = 0). On left panels, the 
dashed and solid curves correspond respectively to the Gaussian limit and our semi-analytic expression (47). On right panels, the solid 
curves correspond to the Gaussian limit. The symbols with errorbars are the measurements. 



n=0). The results are shown in Figure 6 for the differen- 
tial skeleton length (left panels) and for the measured pdf 
(right panels). The errorbars are obtained by the scatter 
over the 100 realizations. They are of the same order for the 
skeleton as for the pdf, although slightly larger for the for- 
mer than for the latter. These small differences will be more 
visible and explained in the next section, which deals with 
non Gaussian cases. 



The case n = requires a larger smoothing window, compared 
to the pixel size, see discussion in Appendix. 10 pixels is a rather 
conservative but safe choice for n = 0. It implies to use images of 
2048 X 2048 pixels, in order to preserve the ratio between the size 
of the smoothing window and the total size of the image, chosen 
for all values of n to be approximately equal to 1/200. 



We see that the departure of /(x,7) from a Gaussian 
(dotted curve on all the panels) is quite weak, even for n = 0, 
which corresponds to a large value of 7 = 0.71. As a result, 
only first order correction is needed, and we find numerically 

that 

Co = 1, C2 = 0.17, (45) 
C2n = for n > 2, (46) 

provides an excellent approximation to /(x,7) in the Gaus- 
sian limit (solid curve on each left panel). Our final expres- 
sion for the normalized differential length of the total local 
skeleton, is thus, for a smoothed Gaussian field: 
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Figure 7. The measured differential length (left panel) of the total local skeleton [eq. (8)] and the pdf of the smoothed field (right panel) 
for distributions with = 1, 2, 3 and 4 degrees of freedom as indicated on the upper right part of each panel. 



3.2 The non Gaussian case 



We now consider a few non Gaussian experiments. The first 
one is the case of a distribution with k degrees of freedom 
(using scale-free Gaussian seeding fields with spectral index 
n = —2). Fig. 7 is similar to Fig. 6, but for distributions 
with k — 1, 2, 3 and 4. The number of realizations, the 
resolution of the maps and the smoothing are the same as 
in Fig. 6. Clearly, the non Gaussian nature of the field is well 
reflected by the skeleton. In fact, and quite surprisingly, its 
differential length again scales very much like the pdf of the 
smoothed density field: there is very little difference between 
left and right panels of Fig. 7, except maybe for the size of 
the errorbars: those are slightly larger for the skeleton than 
for the pdf, as expected. Indeed, the skeleton construction 
relies on estimates of derivatives in a one dimensional subset 
of pixels in the density map, it is therefore more sensitive to 
noise. 



Our second non Gaussian experiment is the Zel'dovich 
map studied in § 2 (Fig. 4). Figure 8 compares the differen- 
tial length of the total local skeleton to the measured pdf. 
Again, for this single realization of a strongly non Gaus- 
sian field, the agreement between both measurements is very 
good, even in the high density tail. Note that the curve for 
the skeleton is slightly more irregular than the one for the 
pdf, as expected. 



Finally, to confirm the validity of the striking results of 
this section, we decided to perform a quite extreme test as 
illustrated by Fig. 9. Taking the field generated on left panel 
of Fig. 1, we increased locally the density contrast by a factor 
1600 along 5 lines 800 pixels long and 1 pixel large, with ran- 
dom positions and random orientations. This map was then 
smoothed with a Gaussian window of radius 30 pixels, as 
shown on upper panel of Fig. 9. The total local skeleton ob- 
tained from this map is displayed in middle panel. The lines 
are clearly visible on the picture as comb like structures. As 
shown in lower panel, the skeleton differential length again 
scales very much like the pdf. 




Figure 8. Comparison of differential length of the total local 
skeleton [points of space verifying eq. (8)] with the pdf of the 
smoothed density field for the Zel'dovich map of Sect. 2, both 
as functions of the normalized density contrast 77 = 6 /a. The 
thick/thin curve corresponds to the skeleton/pdf. For reference, 
the Gaussian limit is also plotted as a dotted curve. 



4 DISCUSSION AND LINKS TO DYNAMICS 

In this paper, we studied some properties of the skeleton of a 
2D random smooth field. This latter is given by an ensemble 
of special field lines connecting saddle points to extrema. It is 
aimed to give accurate account of the network of filaments in 
the field. The skeleton is a nonlocal object, difficult to build 
with a reliable algorithm and to use to perform analytic 
calculations. We thus tried to find a local approximation to 
it, depending on the field and its derivatives. To do that, we 
used two approaches, a mathematically motivated one based 
on a Taylor expansion of the field around local maxima and 
saddle points, and a physically motivated one based on an 
examination of isocontour lines. They both lead to the same 
conclusion. After comparing the resulting local skeleton to 
the real one, we performed statistical analyses of its length as 
a function of density threshold. To simplify the calculations, 
we considered a larger set of curves than the local skeleton 
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Figure 9. The case of a Gaussian field with fines superposed on 
it. 

Top panel: the Gaussian field of left panel of Fig. 1, where the 
density contrast has been increased by a factor 1600 along five 
lines (800 pixels long and one pixel large) randomly located and 
oriented, then smoothed with a Gaussian window of radius 30 
pixels. 

Middle panel: the total local skeleton obtained from the smoothed 
field [points of space verifying eq. (8)]. 

Bottom panel: the measured skeleton length (thick line) compared 
to the pdf (thin line) of the smooth density field, as functions of 
the normalized density contrast rj = 5/(7. The dots correspond to 
the Gaussian limit. Note that due to our procedure, the field has 
very strong kurtosis, but no skewness. Note as well that the cen- 
tral peak is a Gaussian with the variance of the initial Gaussian 
field of right panel of Fig. 1. 



but still containing it, that we called the total local skeleton. 
Having initially in mind to use the total local skeleton as 
a test of non Gaussianity in CMB maps, we compared its 
differential length as function of the density threshold with 
the measured probability distribution function (pdf) of the 
smoothed field. The results of our paper can be summarized 
as follows: 

(i) By definition the real skeleton is the ensemble of pairs 
of field lines departing from saddle points, aligned initially 
with the major axis of local curvature (corresponding to the 
largest eigenvalue of the Hessian) and connecting them to 
local maxima. These field lines are drawn by going along the 
trajectories with the equation of motion dr/dt = Vp until 
convergence to a local maximum. 

(ii) A very good approximation to the real skeleton, the 
local skeleton, is given by points of space where the gradient 
is aligned with the major axis of local curvature and where 
the second component of the local curvature is negative (i.e. 
the smallest eigenvalue of the Hessian is negative). We no- 
ticed however that the local skeleton was shorter than the 
real one, as expected, due to the Lagrangian nature of the 
latter, which can have more than two fields lines converging 
to local maxima, at variance with the former. 

(iii) The total skeleton is given by all the points of space 
where the gradient is aligned with one axis of the curvature. 
Its differential length, as a function of density threshold, is 
seen to scale very closely like the pdf of the smoothed density 
field. This is explicitly demonstrated in the Gaussian case 
with analytic calculations. 

The result mentioned in last point might be discouraging 
for using the total local skeleton as a test of non Gaussian- 
ity in 2D maps, since it does not do better than the much 
simpler pdf. Moreover, since the skeleton depends on local 
derivatives of the density field, it is expected to behave less 
well than the pdf with respect to the noise, although we 
did not investigate that in details, except in part for cosmic 
variance effects. However, our analysis was not exhaustive. 
There might exist some counterexamples where the skeleton 
differential length scales differently from the pdf. Our anal- 
yses relied on isotropic smoothing with a Gaussian window 
and we suspect that this contributes to making the skeleton 
differential length very alike the pdf. More importantly in 
our analyses, we kept the full skeleton length, £(— cxo), as an 
adjustable parameter. In fact, this length does bring signif- 
icant additional pieces of information and should be taken 
into account while comparing models predictions to mea- 
surements. Its analytical calculation in the Gaussian limit, 
although theoretically feasible, is rather cumbersome: we left 
it for future work, having in mind that it can be fairly de- 
termined numerically through appropriate realizations of the 
models. 

Along this paper, we did not address dynamics, al- 
though we used a Zel'dovich map to illustrate our purpose. 
The skeleton is in fact a quite useful tool for the analysis of 
the large scale structure distribution and the understanding 
of its dynamics. Indeed, recall that we initially defined the 
skeleton as the border of void patches (§ 2.1) and that the 
void patches are given by the regions of space containing all 
the points converging to the same local minimum while go- 
ing along the field lines in opposite direction to the gradient. 
On can similarly associate the peak patches to local max- 
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ima (Bond & Myers 1996a). Together with the skeleton, the 
peak patches are the building blocks of the observed large 
scale structures in the Universe. 

Indeed, in the standard approach of hierarchical forma- 
tion of galaxies, the peak patches collapse and merge to- 
gether at successive times. The merging of peak patches can 
also be seen as the collapse of larger peak patches consti- 
tuted by their progenitors. Within a multi-scale approach, 
these latter can be obtained by smoothing the field at in- 
creasing scales, each smoothing scale corresponding to a dif- 
ferent collapse time, as illustrated by Fig. 10. This peak 
patch approach was in fact used extensively by Bond & 
Myers ( 1996a, b,c) to produce a simplified but quite accu- 
rate description of large structure dynamics, including the 
evolution of non linear objects from galaxies to clusters of 
galaxies, their merging history as well as their large scale 
motions. Note that peak patches are rather compact, and 
thus can be well approximated by ellipsoids, lending cre- 
dence to the Bond & Myers approach. This also stems from 
a Taylor expansion of the field around local maxima. 

This line of thought can be followed further. Indeed, 
the local maxima are by definition located on the skeleton 
along which the matter fiows: merging of earlier collapsed 
patches will take place at the nodes of the skeleton (see also 
caption of Fig. 10). This description is well known and un- 
derstood through the adhesion approximation (e.g., Kofman 
& Shandarin 1988; Kofman, Pogosyan & Shandarin 1990). 

Note however that the skeleton itself has its own dy- 
namics, as illustrated by Fig. 11: filaments composing it can 
move and be distorted due to large scale fiows, but can also 
merge together. For the particular example considered here. 
Fig. 5 shows that the total length of the skeleton is approxi- 
mately conserved. It is slightly shorter in the Zel'dovich case 
compared to initial conditions as a result of the competition 
between merging and stretching. 

Since the matter tends to fiow along the lines of the 
skeleton, these latter represent a useful reference frame to 
study internal structure and dynamics of filaments in the 
Universe, with well prescribed procedure to define them, 
given a typical scale length. In practice, this latter should be 
larger than the size of clusters of galaxies, and smaller than 
the size of super-clusters, i.e. of order of a few Mpc. In that 
case, one knows that the skeleton refiects rather well the ini- 
tial primordial one (Bond, Kofman & Pogosyan 1996). Note 
again that the skeleton length can be measured in a cosmo- 
logical volume, and compared to theoretical predictions. 

Since we are in the 2D case, the discussion concerning 
the dynamics remained at the qualitative level. More quan- 
titative analyses in 3D are left for future work. 
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Figure 10. The skeleton and the peak patches for the Gaussian field of Fig. 1, at two different smoothing radii i. 

The left panel corresponds to -^i =25 pixels, as previously. The right panel corresponds to -^2 = 25 a/2 ~ 36 pixels. The skeleton is 
represented by the blue lines, and the peak patches borders (the dual skeleton) by the golden ones. Each red point corresponds to a 
local maximum. Left and right panels can be seen as a Lagrangian view of the system at two different times, describing the merging of 
structures. 

Clearly, the peak patches of the right panel can be seen as mergers of peak patches of the left panel, even if mergers can occur differently 
according to the place of interest: some peak patches survive, i.e. do not merge with others, some of them experience merging with one 
or more neighbors. As a result, the skeleton of right panel is approximately made of a connected subset of lines composing the skeleton 
of left panel. 

Note that on the lower left corner of left panel, there is a peak patch containing no local maximum. This is clearly an artifact from our 
numerical approach, which did not detect it. This is not surprising since we noticed earlier (Fig. 2) that something was wrong with the 
connectivity of the skeleton between critical points in this location. 




Figure 11. Dynamical evolution of the skeleton for the smooth Gaussian field of Fig. 1. 

Left panel: the final skeleton obtained from the Zel'dovich mover (i.e. as explained in Fig. 4, golden lines) superposed to the initial one 
(blue lines). 

Right panel: same as left panel but for over-dense regions. 

One can clearly match the details in the initial pattern of filaments to the final one, except for mergers as in e.g. the top right corner 
and the left side of the panel. Note as well the expansion of the central under-dense region. 

tensively pixelization and finite volume effects by generating 
maps at various resolutions and smoothed at various scales. 
We tried as well various schemes described below for com- 
puting derivatives and interpolating the field, and we run 
many different realizations of the same power-spectrum. Fi- 
nite volume effects are more important for smaller values of 
the spectral index, n, while pixelization effects, on the con- 
trary, increase with n. In principle if the smoothing length 
I is very large compared to the pixel size and very small 



APPENDIX A: NUMERICAL APPROACH 

The numerical calculations were performed with a dedicated 
FORTRAN 90 package: FasToChtt (Fast Topological Chase). 

While confronting the measurements to the theoretical 
predictions of Sect. 3.1 for the Gaussian case, we tested ex- 

M Available on request from the authors. 



© 0000 RAS, MNRAS 000, 000-000 



Skeleton as a probe of the cosmic web 15 



compared to the map size, L, both these effects should be 
neghgible, as we found for n = 0, — 1 and —2. Practicahy, 10 
pixels ^ £ ^ L/20 is generally a safe choice, but in fact the 
range of available scales depends on the statistics considered 
as discussed more in detail below an on details of the field 
properties, in particular its power-spectrum shape. We also 
tested anisotropy effects by generating the initial random 
field with or without Hanning filtering, the latter insuring 
isotropy at small scales (e.g. Berstchinger 2001) and found 
that they were insignificant. 

Before entering into the details of the skeleton con- 
struction, we first detail the issue of computing reliably 
the successive derivatives of the field. We examined both 
Fourier methods and the simplest finite difference schemes. 
Throughout this paper, we used the Fourier method but the 
simple finite difference method is both faster and easier to 
use when the coverage is more intricate such as in galaxy 
catalogs or in CMB experiments. We also investigated pix- 
elization effects and found that they are neghgible provided 
that the smoothing window is large enough, i.e. a few pixels 
radius, typically £ ^ 3. However this of course depends on 
the type of unsmoothed map considered: if there is a lot of 
small scale power in the map, it is necessary to smooth it 
more to have reliable estimates of the derivatives. 

We also tested bicubic interpolation (e.g. Press et al. 
1992) which has the advantage of warranting divergence- 
free gradient at all positions within a pixel. Practically, this 
means that given the field values at the pixel centers, as 
well as pre-computations of the gradient and of the diagonal 
terms of the Hessian (with either Fourier or finite difference) , 
any quantity can be computed self-consistently at any loca- 
tion within the pixels, including the off-diagonal terms of 
the Hessian. This property is in principle particularly criti- 
cal when building the real skeleton. In practice, however, the 
improvements brought by the bicubic interpolation were in- 
significant as compared to the simpler and faster bilinear 
interpolation, which was finally used for all the calculations. 

Drawing the local skeleton is simple if one sees that 
the equation S = [eq. (8)] corresponds to the zero isocon- 
tour of the field S. This can be performed with a standard 
method as we now explain. Given a square of four neigh- 
boring pixels and the corresponding values of <S, we first de- 
termine whether the linearly interpolated field cancels along 
two sides of the square. If this happens, this means that 
the isocontour curve crosses the square. We locally approxi- 
mate this curve by a segment which extremities are located 
on the edges of the square. The coordinates of the segment 
ends are easily found using dual interpolation. The length 
of the skeleton is found by adding all the individual segment 
lengths. There are however particular cases where two pieces 
of isocontour can intersect at the same point (e.g., at a crit- 
ical point) or become very close to each other. This can pro- 
duce configurations where the field cancels on 4 edges of the 
square. In that case, we cannot compute reliably the length 
of the isocontour within the square, but the relative con- 
tribution of these configurations is increasingly small with 
the smoothing length. To make these contributions negligi- 
ble, and to insure as well that the effect of approximating 
the local skeleton locally by straight lines is negligible, the 
smoothing length should be of order a few pixels size, typi- 
cally £ ^ 5 for spectral index n ^ — 1. However, the smooth- 
ing radius has to remain small compared to the map size. 



in order to avoid finite volume effects, typically i ^ L/20, 
where L is the map size. 

Drawing the real skeleton is rather difficult, at least we 
did not find yet a highly reliable algorithm. As explained 
in Sect. 2.1 the real skeleton is drawn by going along the 
trajectory with the following motion equation 

|.v = Vp, (Al) 

starting from the saddle points and with initial velocity par- 
allel to the major axis of the local curvature (the dual skele- 
ton is obtained similarly by making the operation p —p). 
This equation of motion is solved numerically with a semi- 
implicit scheme which guaranties that the gradient does not 
change sign along the direction of motion. This allows as 
a consequence a correct calculation of the skeleton length 
since backwards motion along the trajectory is prevented. 

More explicitly, if r^, = Vp* and dti are respectively 
the position, velocity and timestep at step z, the quanti- 
ties at next step are computed as follows using an itera- 
tive procedure. The gradient Vp is evaluated at the po- 
sition fi+i = n -\- Vi-^idti-^i from which we deduce v^+i 
and s = Vj+i.Vi. As a first guess we take dii-^i = dU and 
Vi+i = Vi. As long as s < we divide the timestep dtiJ^i 
by two and recompute fz+i. Once s > 0, the quantities at 
time step i + 1 are set to dU^i = c/t^+i, v^+i = v^+i and 
1*2+1 = + Vi+idti+i. Note in addition that the time step is 
always chosen such that the displacement between two steps 
is at most of the order of a pixel. 

To ensure a finite number of time steps, the motion is 
stopped once we reached the vicinity (typically a fraction 
of pixel size) of a critical point and appear to converge to- 
wards it. In principle this stopping point should always be 
a maximum. However, in practice we found that very rarely 
we could converge to a saddle point. To avoid here again 
an infinite number of iterations, we stop the calculation of 
the trajectory. Indeed, this situation happens when a field 
line (say A) is very close to one of the two unstable field 
lines converging to the saddle point (see Fig. 3). After get- 
ting close to the saddle point, the field line A should turn 
by approximately 90 degrees to follow closely one of the two 
stable field lines (say B) of the saddle point. The field line 
B will be drawn anyway starting from this saddle point. So 
visually, we should not miss anything by not drawing the 
end of the field line A. However, this might be a problem for 
estimating the skeleton length since we are missing parts of 
it. 

Eventually, the critical points are determined as inter- 
sections of the contour lines dp/dri = and dp/dr2 = 
with the same method as used for the local skeleton. This 
detection method can also fail (e.g. upper left panel of Fig. 2 
and 10), and some critical points can be missing, which has 
dramatic consequences for the building of the real skeleton. 

To avoid critical situations as described above, where 
our algorithm for drawing the real skeleton fails, we need 
to smooth the fields significantly, typically with a radius at 
least of order of 10 pixels but this depends strongly on the 
nature of the field. 
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